
# A3c_price_index_full_adj
#==============================================================================

# Description: This code computes the difference in the gains from trade between
# poor and rich on average across sectors in the full model

rm(list = ls())
options(warn=-1)

library('data.table')
library('stargazer')
library('ggplot2')
library("ggrepel")  


# Import product-specific price indexes
#===============================================================================

i_k = 1
data_name <- sprintf("price_index/product_specific/diff_full_adj_table_%s_k.csv", i_k)
data = read.table(data_name, header =TRUE, sep=",")
data <- as.data.table(data)

# Create index that is 1 if |beta| < 0.6 or if the Home country's mc could not be inferred (demand not sufficiently elastic)
data <- data[ , index_mc := (any(abs(beta_price_80) < 0.6) | any(abs(beta_price_20) < 0.6) | any(Declarant == Partner & mc_flag == -2)), by = c('Year', 'Quarter', 'Declarant')]

# Keep 1 observation for each market 
data <- data[ , .SD[1], by = c('Year', 'Quarter', 'Declarant')]                               
data <- data[, c('Year', 'Quarter', 'Declarant', 'P_true_sorted_20','P_true_sorted_80', 'P_counter_sorted_20','P_counter_sorted_80', 'beta_price_20', 'beta_price_80', 'index_mc'), with = FALSE]

data <- data[ , product := i_k]   

I_K = 3426
for (i_k in 2:I_K) {
  
  print(i_k)  
  
  # Import and Clean Data  
  data_name <- sprintf("price_index/product_specific/diff_full_adj_table_%s_k.csv", i_k)
  data_temp = read.table(data_name, header =TRUE, sep=",")
  data_temp <- as.data.table(data_temp)
  
  # Create index that is 1 if |beta| < 0.6 or if the Home country's mc could not be inferred (demand not sufficiently elastic)
  data_temp <- data_temp[ , index_mc := (any(abs(beta_price_80) < 0.6) | any(abs(beta_price_20) < 0.6) | any(Declarant == Partner & mc_flag == -2)), by = c('Year', 'Quarter', 'Declarant')]
  
  # Keep 1 observation for each market 
  data_temp$q_new = NULL
  data_temp <- data_temp[ , .SD[1], by = c('Year', 'Quarter', 'Declarant')]      # Keep 1 obs for each market 
  data_temp <- data_temp[, c('Year', 'Quarter', 'Declarant', 'P_true_sorted_20','P_true_sorted_80', 'P_counter_sorted_20', 'P_counter_sorted_80', 'beta_price_20', 'beta_price_80', 'index_mc'), with = FALSE]
  
  data_temp <- data_temp[ , product := i_k]  
  
  data = rbind(data, data_temp)
  
}  


save.image('data_PI_full_adj.RData')
#load('data_PI_full_adj.RData')

#===============================================================================


# Format
#===============================================================================

data[ , u_beta_true_20 := log(P_true_sorted_20) - 1]                                     
data[ , u_beta_true_80 := log(P_true_sorted_80) - 1]                                    
data[ , u_beta_counter_20 := log(P_counter_sorted_20) - 1]
data[ , u_beta_counter_80 := log(P_counter_sorted_80) - 1]

# Create difference between poor and rich
#-------------------------------------------------------------------------------

data[ , diff := P_counter_sorted_20 / P_true_sorted_20 - P_counter_sorted_80 / P_true_sorted_80]

# Remove extreme observations
#-------------------------------------------------------------------------------

data = data[diff < 1 & diff > -1]
data = data[(is.na(u_beta_counter_80) | is.na(u_beta_counter_20) | is.na(u_beta_true_80) | is.na(u_beta_true_20)) == F]
data = data[(abs(u_beta_counter_80) < 1e-13 | abs(u_beta_counter_20) < 1e-13 | abs(u_beta_true_80) < 1e-13 | abs(u_beta_true_20) < 1e-13) == F]
data = data[index_mc == F]                                                      # Remove cases in which demand is not sufficiently elastic 

data[ , v_true_20 := beta_price_20 * log(P_true_sorted_20 / (exp(1) * gamma(1 - (1/beta_price_20))))]
data[ , v_counter_20 := beta_price_20 * log(P_counter_sorted_20 / (exp(1) * gamma(1 - (1/beta_price_20))))]
data[ , v_true_80 := beta_price_80 * log(P_true_sorted_80 / (exp(1) * gamma(1 - (1/beta_price_80))))]
data[ , v_counter_80 := beta_price_80 * log(P_counter_sorted_80 / (exp(1) * gamma(1 - (1/beta_price_80))))]

data = data[(v_true_20 == Inf | v_counter_20 == Inf | v_true_80 == Inf | v_counter_80 == Inf) == F]
data = data[(v_true_20 == -Inf | v_counter_20 == -Inf | v_true_80 == -Inf | v_counter_80 == -Inf) == F]

#===============================================================================


# Run Regression, save, and plot
#===============================================================================

reg_full_adj <- lm(diff ~ 0 + as.factor(Declarant), data = data)
print(summary(reg_full_adj))

declarant <- as.matrix(reg_full_adj$coefficients)
declarant <- gsub("as.factor\\(Declarant\\)","",rownames(declarant))
declarant <- as.numeric(declarant)
income_data <- fread(file = "gdp_per_capita_2007.csv")

# Save Output

output <- merge(income_data, data.table(declarant, reg_full_adj$coefficients), by = "declarant")
colnames(output) <- c(colnames(output)[1:4], "FE_declarant")
write.csv(file = "price_index/results/diff_full_adj.csv", output)


# Plot
#-------------------------------------------------------------------------------

rm(list = ls())
output <- fread('price_index/results/diff_full_adj.csv')

iso3 <- fread('ISO3_declarant_codes.csv')

output[ , gdp_log := log(gdp_per_capita_declarant)]
setkey(output, declarant)
setkey(iso3, declarant)

output <- output[iso3]
output <- output[is.na(country) == F]
output = output[declarant != 18] 
output = output[declarant != 600] 
output = output[declarant != 46] 


reg <- lm(FE_declarant ~ gdp_log, data = output)
summary(reg)
cor(output$FE_declarant, output$gdp_log)

plot_diff_full_adj <- ggplot(output, aes(x=gdp_log, y=FE_declarant)) +
  geom_text_repel(label=output$ISO3, size = 5) +
  theme(text = element_text(size=17)) +
  geom_smooth(method=lm) +
  labs(x = "GDP per capita (in logs)", y = expression(paste(Delta, "Welfare"["Poor"]," - ", Delta, "Welfare"["Rich"]))) 

ggsave("price_index/results/plot_full_adj_R.png", plot = plot_diff_full_adj)

